One-dimensional gas of hard needles 
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We study a one dimensional gas of needle-like objects as a testing ground for a formalism that 
relates the thermodynamic properties of "hard" potentials to the probabilities for contacts between 
particles. Specifically, we use Monte Carlo methods to calculate the pressure and elasticity coefficient 
of the hard-needle gas as a function of its density. The results are then compared to the same 
quantities obtained analytically from a transfer matrix approach. 
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I. INTRODUCTION 

Due to the relative simplicity of the derivation of their 
thermodynamic properties, classical one-dimensional 
(ID) systems are frequently employed as a test bed of 
theory and methods for collective behavior in higher di- 
mensional systems. For instance, the collection of "hard 
spheres" on a line, sometimes referred to as the Tonks 
gas [l| , has served as an initial step in the study of two- 
and three-dimensional systems of hard disks/spheres. 
There is indeed a general method for exact analysis of 
a gas of point particles interacting in ID via potentials 
that depend only on near- neighbor separations [2]. Here, 
we employ such methods to study a gas of hard needle- 
shaped objects. Our object is to compare the analytical 
results with those obtained by (Monte Carlo implemen- 
tation) of a formalism that relates thermodynamic prop- 
erties (specifically the pressure and elasticity coefficient) 
of the gas to the probabilities of contact amongst the 
particles. 

"Hard" potentials, which are either zero or oo, help to 
illuminate the geometrical/entropic features of a thermo- 
dynamic system. Since there is no energy scale arising 
from such potentials, the temperature T appears only 
as a multiplicative factor in the free energy and various 
other thermodynamic quantities, such as pressure and 
elastic coefficients. Thus the state of the system becomes 
independent of T, and only depends upon such features 
as density. The clarity of the geometrical perspective, 
combined with the simplicity of numerical simulations, 
has lead to extensive studies of such systems. In fact, 
simulations with hard potentials date back to the origins 
of the Metropolis Monte Carlo (MC) method 0, and 
have flourished in the decades that have followed (see 
Ref. [4] and references therein). A typical example of 
non-trivial behavior is the entropically driven first order 
phase transition from a liquid to a solid phase [5[ . 

Alignments of non-spherically-symmetric molecules 
lead to a diversity of phases in liquid crystals [|]. For 
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FIG. 1: (Color online) Needle-shaped particles of length 2£ 
and vanishing thickness are free to rotate in two dimensions, 
with their centers moving along a line. Particle (needle) i is 
characterized by its translational position x% and orientation 
angle fa. 



example, in the nematic phase the molecules have no po- 
sitional order (like a liquid), while their orientations are 
aligned to a specific direction. From the early stages re- 
search into liquid crystal it was realized that the entropic 
part of the free energy related to non-spherical shapes of 
the molecules, can by itself explain many of the proper- 
ties of such systems [7| . Not surprisingly, hard potentials 
were frequently invoked, and even such simplifications as 
infinitely thin disks @ or rods @ have provided valu- 
able insights regarding liquid crystals. The interplay be- 
tween the rotational and translational degrees of freedom 
in molecular solids [10] leads to elastic properties that are 
coupled to orient at ional order. How does one compute 
the elastic response of such systems from first principles? 

Recently a formalism enabling direct calculation of 
elastic properties and stresses of a system of hard 
non-spherically-symmetric objects was developed [ll| 
by extending a previously known formalism for hard 
spheres [12j. Not surprisingly, given that the elastic re- 
sponse in two and higher dimensions depends on a rank 
four tensor, the resulting expressions contain a large 
number of terms. Typical terms correspond to a vari- 
ety of possible contacts between particles and numerous 
components of the separations between them (see, e.g., 
Eq. 23 in Ref. [11]). Since these expressions are ob- 
tained after numerous mathematical transformations, it 
is advisable to subject them to independent tests. In- 
deed they have been shown to reduce to the known re- 
sults for isotropic objects, but up to now there had been 
no comparison with exact results for non-spherical par- 
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tides. Here, we consider the statistical mechanics of a 
ID system of hard needle-like particles rotating in two 
dimensions with their centers affixed to a ID line, as de- 
picted in Fig.[H The needles are not allowed to intersect, 
and thus act as "hard" potentials. This model is a par- 
ticular case of a group models considered by Lebowitz et 
al. [l3[ with anisotropic objects in one dimension. From 
the perspective of complexity, such systems are a slight 
generalization of the Tonks gas, yet they provide non- 
trivial insights into the interplay of rotational and trans- 
lational degrees of freedom. The model can be solved ex- 
actly, and thermodynamic properties, such as elastic co- 
efficients, can be calculated. We compare the values ob- 
tained analytically by the transfer matrix method, with 
those from MC simulations using the expressions from 
Ref. [11] adapted to the ID case. 

The paper is organized as follows: The model of hard 
needles is introduced in Sec. [Ill an d we demonstrate how 
the relative orientations of neighbors leads to an effective 
hard-potential as a function of their separation. Sec. HIJ 
is devoted to reviewing how elastic properties of a system 
can be characterized, and the expressions for computing 
elastic coefficients in ID are presented. The numerical 
difficulties associated with evaluation of various quanti- 
ties by MC simulations are also described. In Sec. HVlwe 
present the transfer matrix method for solution of the 
model. Details of the MC simulation are presented in 
Sec. El along with comparisons to results obtained by 
transfer matrix method. Discussion and additional fea- 
tures of the model are presented in Sec. I VII 

II. MODEL 

Figure [1] depicts a configuration of our model, consist- 
ing of needles of length 21, with their center positions re- 
stricted to move on a ID line. Needle i is characterized by 
its position Xi, and orientation fa measured with respect 
to the normal to the line. Since orientations differing by 
7r are indistinguishable, we restrict — tt/2 < fa < tt/2. 
(Such entities, called directors, frequently appear in the 
description of liquid crystals [6].) As i is the only mi- 
croscopic length scale in the problem, it can be used to 
construct dimensionless parameters. In particular, the 
mean distance between particles a is made dimensionless 
by considering a/ 1, while the density n can be replaced 
by ni. 

The needles are not allowed to intersect but do not 
interact otherwise. Since the particles cannot cross each 
other, we number them (left-to-right) along the ID line, 
and require that this order is unchanged, i.e. Xi-i < 
Xi. (This convention simplifies the enumeration of the 
possible contacts between particles.) Thus the distance 
of closest approach between adjacent needles is a function 
of their orientations, given id, with the dimensionless 
function 

di-l,i{(pi-l,(pi) = ; 7- r 7— -tt . (1) 

max[cos(0i_ij, cos(0ijj 




FIG. 2: Gray-level representation of the function di^ifaifa) 
in Eq. (pQ) for the dependence of minimal distance between 
needles on their orientations. The black diagonal corresponds 
to d = for parallel needles, while white corresponds to d = 2 
for needles along the line. 



This function is depicted in Fig. [21 and varies between 
zero (when the needles are parallel, fa = fa-i) to 2 
(with the needles lying on the line, fa = —fa-i —> k/2). 
Note that value of d is poorly-defined at the points 
(±7r/2, ±7r/2), and depends on the limiting procedure. 
Analytic computations would have been considerably 
simplified if d was only a function of the difference in 
orientation, but this is not the case because of the de- 
nominator in Eq. (pQ). 

We consider a collection of N needles, either in an en- 
semble of fixed length L (for MC simulations), or fixed ex- 
ternal pressure (force) p (for transfer matrix studies). It 
is convenient to impose the boundary conditions through 
the definition of the minimal distance. For the MC 
simulations, we introduce fictitious particles i = an 
i = N + 1. Periodic boundary conditions on a line of 
length L are implemented by requiring xn+i = x\ + L 
and (j)N+i = 01, while xo = xn — L and fa = c/)n- This 
extends the validity of Eq. (pQ) to i = 1 and i = N + 1, 
enabling the treatment of all particles on equal footing. 
In the fixed pressure ensemble, which is used in trans- 
fer matrix calculations, the orientation of the first (last) 
particle i = 1 (i = N) is restricted only by its neighbor 
from the right (left), i.e. i = 2 (i = N — 1). The position 
of the both end particles is arbitrary in this ensemble, 
with x\ > and xn = L (which is also a variable in this 
ensemble). In this case, we set c?o,i = ^at,at+i = 0. 

As explained above, adjacent particles interact via the 
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hard potential 

Vi- lti = V(x t 



if X{ Xi—\ ^> ^di — X^ii 

otherwise. 



(2) 



As such a potential does not have an energy scale, the 
temperature T will appear only as a prefactor in the ther- 
modynamic quantities. In particular, the Helmholtz free 
energy F which is an extensive quantity with units of 
energy will have a form F = NkBTh(n£), where k B is 
the Boltzmann constant. In ID, the pressure p, and the 
elastic coefficient C have units of force and can be made 
dimensionless by considering / = (3pi and (3C£, where 
(3 = l/k B T. The Gibbs free energy G = Nk B Tg(f) 
depends only on the dimensionless pressure /. 



elastic moduli are related to thermal averages of products 
of various inter-particle forces and separations. This for- 
malism was extended to hard potentials in Refs. pH.H3|. 
Since for hard potentials the forces vanish except when 
the particles touch, the results depend on various contact 
probabilities. In two and three dimensions the stress and 
the elastic constants are tensors, and the expressions in- 
volve averages over a variety of components. These re- 
sults simplify in ID, and in particular, the expression for 
stress [Eq. (22) in Ref. [ll|] can be considerably simpli- 
fied: If we denote the separation between two adjacent 
needles by Si = Xi+i — Xi, they are in contact if the ar- 
gument of Ai = S[si — idij+ifyi, <^i+i)] vanishes. The 
dimensionless pressure then / becomes 



/ = nt 1 + 



N ^ 



siAi) 



(6) 



III. ELASTICITY OF ID 

Shape and size deformations of objects are usually de- 
scribed by the strain tensor [141 ]. In ID this reduces to 
a scalar quantity n which simply relates the distorted 
size of the system V to its original size L via L' 2 = 
L 2 (I + 277) uM- (While this definition is slightly awkward 
in ID, the use of squared distances between points is con- 
venient because in higher dimensions it clearly separates 
trivial changes in geometry caused by rotations, and real 
deformations.) In ID, for small 77, the Helmholtz free 
energy can be expanded as 



(3) 



where C is the elastic coefficient of the body. (Note, 
that the free energy on the left hand side (l.h.s.) of the 
equation is divided by the undistorted size of the system.) 
Consequently, p and C can be calculated from the first 
and second derivatives of F with respect to n at fixed T. 

While the elastic properties are more naturally ob- 
tained from the Helmholtz free energy, we will also use 
the Gibbs free energy G = F + pL, in which the pres- 
sure is the (imposed) variable [16]. The system size L, 
or the mean inter-particle distance is then obtained from 
a = d(G/N)/dp\ T , or in terms of dimensionless variables 



Similarly, C 



a _ difiG/N) 
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Squire et al. [17[ developed a formalism for a direct cal- 
culation of elastic parameters from the correlation func- 
tions of particles. In this approach stress (pressure) and 



The first term in this expression is simply the pressure 
of the ideal gas, while the second term can be easily rec- 
ognized as the mean value of the product of the inter- 
particle separation and force, as appears in the virial the- 
orem [18]. To evaluate Eq. ([6j), we need the probability 
that two particles (i and z + 1) with specified orientations 
(cj)i and touch each other. 

Similarly, the elastic coefficient C can be expressed as 



(3C£ = ni 



-—^(siSjAiAj) 

i<j 



— Y(( 

2N ^ u 



(7) 



5?+5? +1 )AiA<+i) 



The last two sums in the right hand side (r.h.s.) of 
Eq. (0) involve averages of products of As, i.e. they re- 
quire knowledge of the joint probability density of two 
simultaneous contacts. The last sum involves cases when 
three particles z, i+1 and i + 2 touch each other, while the 
preceding sum depends also on cases when two indepen- 
dent pairs are in contact, i.e. particle i touches z + 1, and 
a different particle j (> i + 1) touches j + 1. The l.h.s. of 
Eq. (|ZJ) is an intensive quantity, while the third and the 
fourth terms on its r.h.s. contain 0(N 2 ) terms. How- 
ever, most of the terms appearing in these two sums can 



be grouped in pairs (siAi)(sjAj) 



(siSjAiAj) 



which 



decay to zero when the distance between the pairs of 
particles exceeds the correlation length. All the averages 
appearing in Eqs. (|6|) and can be calculated in MC 
simulations. 



IV. TRANSFER MATRIX APPROACH 

The partition functions of ID models with short range 
interactions can be found analytically using a transfer 
matrix method [l6|, [l9|, [2(j. It convenient to consider 
the isobaric ensemble with fixed external pressure (force) 



p, such that (the configurational part of) the partition 
function is given by 

r N N 

Z G = ]J d ^U d ^~ 0i: - lV '- 1 -~ 0PXN - ( § ) 

1=1 1=1 

Since xn = ^ s ii we can change variables and perform 
integrations over the separations between adjacent par- 
ticles. For the hard potential given by Eq. ([2]), this leads 
to 

= / nf*^-!^/;^!,*)), (9) 

i=i 

where 

A- M (/;^-i,^) ^e-^ 1 - 1 -^- 1 -**). (10) 

(According to the definition of all D are identical, ex- 
cept for Do,i = Dn,n+i = 1 at the boundaries, as ex- 
plained in the Sec. fTTl) 

The expression in Eq. (pQ) is too complicated for the 
integrals in Eq. ([9]) to be performed analytically. Nev- 
ertheless, multiple integrals of this kind can be easily 
performed numerically to any desired accuracy. We can 
subdivide the range of the angular integration into M 
equal segments, by setting <pk = — 7r/2 + -^/c, with 
k = 0, 1, • • • , M — 1. This replaces the function D by 
an M x M matrix, and the integrals in Eq. ([9]) are re- 
placed by matrix products. The partition function then 
becomes 

Z G = (f3p)- N (7r/M) N vD N - 1 v, (11) 

where v is a column vector with all of its elements equal 
to 1. Repeated multiplications can be performed numer- 
ically, first multiplying D by itself, then multiplying the 
resulting matrix by itself, etc.. After a total of K such 
iterations we arrive at D N , with N = 2 K + 1. The expo- 
nential dependence on K allows us to achieve very large 
values of TV, in practice we used K = 20 in our simula- 
tions. For moderate pressures, the discretization of the 
angle (j) has little influence on the result once M exceeds 
10, and we report results for M = 512. (It should be 
noted, that the same results can also be obtained by nu- 
merically finding the largest eigenvalue of D. However, 
in our case this alternative provides no numerical advan- 
tage.) From the numerical value of Zg, we then obtain 
the Gibbs free energy 

Figure [3] depicts the scaled Gibbs free energy calcu- 
lated by this numerical procedure. For non-interacting 
needles the partition function is Zq = (7r//3p)~ N , and 
the corresponding /3Go/N = \n(/3p/7r) is indicated by the 
lower line in the figure. (Both curves exclude the trivial 
contribution due to kinetic energy.) 




FIG. 3: (Color online) The upper curve depicts the Gibbs 
free energy per particle (made dimensionless by multiplying 
by 0) as a function of the dimensionless pressure / = f3p£. 
The lower curve shows, for comparison, the same quantity for 
non-interacting needles; the curves begin to separate when / 
is larger than about one. 



V. SIMULATIONS AND RESULTS 

A Monte Carlo procedure was used to evaluate the 
pressure and elastic coefficient of the system of hard nee- 
dles. We simulated TV = 128 particles with periodic 
boundary conditions. Correlations between the needles 
for small and moderate densities do not persist past a 
few neighbors, and our choice of N caused no percepti- 
ble finite size effects. An elementary MC move consists of 
randomly choosing a needle, randomly deciding whether 
to displace or to rotate it, and attempting to perform 
such a move. The move is accepted if in the new posi- 
tion, or with new orientation, the needle does not overlap 
with neighboring needles. The particles are sequentially 
ordered, and a position change is rejected if it changes 
this ordering. The attempted moves are uniformly dis- 
tributed over an interval, whose width is chosen to be as 
large as possible, while maintaining an acceptance rate 
larger than 50%. The varying size of the interval implies 
a diffusion constant for each particle that decreases with 
increasing density. A single MC time unit consists of 
2N attempts to move or rotate particles. The relaxation 
time of the system is proportional to L 2 , and inversely 
proportional to the diffusion constant and elastic coef- 
ficient. (The latter increases with increasing density.) 
Our choice of elementary step ensured that, within the 
examined range of densities, the relaxation time was ap- 
proximately constant and remained of order N 2 . This 
was verified by directly measuring several autocorrela- 
tion functions. For every density n the simulation time 
was 5 x 10 5 A/" 2 . Such long times are required to ensure 
high accuracy of measured contact probabilities, as ex- 
plained below. 

The presence of the Dirac ^-function in the definition 
of Ai necessitates delicate handling. Both Eq. (|6]) and 
Eq. require measuring separation Si between the adja- 
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FIG. 4: (Color online) The mean distance between particles 
a (in units of fl, as a function of dimensionless force /, for an 
ideal gas (lower continuous line) as compared with the virial 
expansion truncated at the second term (upper continuous 
line), and with the exact transfer matrix result (dashed line). 
The dashed-dotted line is the transfer matrix result for the 
elasticity coefficient C. Solid circles represent the relation 
between a and / from the MC simulations. Solid squares 
represent the MC results for the elastic coefficient. 



cent needles at the moment of contact. Such events have 
zero probability, and the formulas really involve probabil- 
ity densities. The latter can be evaluated by examining 
the probability that the two needles are within ei, and 
divide the result by e\. Of course, the number of such 
near collisions decreases with decreasing ei, and the sta- 
tistical error increases. The situation is even worse for 
terms of type (s^SjA^Aj), where two simultaneous con- 
tacts are supposed to appear. One may define two near 
collision events by considering intervals of size e\ and 
62- The opposing requirements of having — > (for 
accurate calculation of probability densities), and large 
€i (to ensure statistical accuracy) can be partially recon- 
ciled by considering each argument of a (5-function being 
in the range [me, (m + l)e), with m = 0, 1, • • • , M. We 
used M = 10, and e = 0.002 (0.01) for high (low) par- 
ticle density simulations. With 11 data points for single 
contact terms, and ll 2 points for two-contact terms, we 
could view the results as a function of one variable mi, 
or two variables mi and m2, and extrapolate the results 
to the "exact contact" limit. The accuracy and prac- 
ticality of such a procedure has been demonstrated in 
Ref. [12]. The total simulation time was determined by 
requirement of having sufficient number of terms in each 
"bin" of the statistical procedure explained above. The 
total simulation time was dictated by the need to have 
a very accurate estimate of the fourth term on the r.h.s. 
of Eq. ©. 

Since the MC simulation is performed in the ensem- 
ble of fixed length, the density or mean inter-particle 
distance a are given, while the dimensionless pressure 
/ and the dimensionless elastic modulus /3£C are calcu- 
lated. The full circles in Fig. 2] depict the calculated 



dependence of / (horizontal axis) on a (vertical axis). 
The error bars on / are negligible, since Eq. (j6| includes 
only single pair contacts, and the large statistics as well 
as small e\ ensures very high accuracy. This result is 
compared with the relation between / and a obtained 
from the Gibbs free energy via Eq. (jlj) by taking the nu- 
merical derivative of G calculated by the transfer matrix 
method. The latter is depicted by the dashed line. Ex- 
cellent agreement is obtained between the results from 
these two methods. 

The solid squares in Fig.JUdepict the MC results for the 
dimensionless elastic coefficient /3£C (vertical axis), as a 
function of the dimensionless pressure / (horizontal axis). 
Since in the MC procedure / is itself a computed quan- 
tity, there are now also horizontal error bars, which are 
negligible as explained in the previous paragraph. The 
accuracy of C, however, is much lower and depends on 
both statistical errors, and systematic errors from ex- 
trapolation to the true contact probability densities. We 
chose the values of and the simulation time in such a 
way that both errors were of the same order. We esti- 
mate that the vertical error bars are approximately the 
size of the symbol for the leftmost point, and decrease to 
half the symbol size for the rightmost point. The dashed- 
dotted line depicts the same relation obtained from G by 
using Eq. ([5]) and the transfer matrix calculations. The 
results from both approaches coincide within the esti- 
mated errors. 



VI. DISCUSSION 

The good agreement between the results from MC sim- 
ulations, based on contact probabilities, and those from 
the transfer matrix method, support the validity of the 
expressions reproduced in Eqs. (|6]), ([7j) for the pressure 
and elastic moduli of hard potentials. While limited to 
ID, this is the first direct comparison between formulae 
derived in Ref. [ll|, and an exact alternative approach. 

We conclude by pointing out an interesting feature of 
the hard needle system: Both at small and large den- 
sities, the pressure and elastic coefficient are related by 
the simple expression C = 2p, while the behavior at in- 
termediate densities is more complicated. For low pres- 
sure (density) the system behaves as an ideal gas with 
a = (/?£>) -1 , and substituting in Eq. ([5]) immediately 
yields C = 2p in this limit. Interestingly, as discussed 
by Lebowitz et al. [13], the relation between density and 
pressure also simplifies at very high pressure (density). 
In this limit the angular integrations are themselves con- 
strained by pressure, and a Gaussian approximation leads 
to additional powers of /3p in the Gibbs partition func- 
tion. This in turn leads to a density n = a -1 = 2(3p, i.e. 
the same functional dependence as an ideal gas but with 
a factor of 2. Inserting this limiting behavior into Eq. ([5]) 
again leads to C = 2p as in the ideal gas limit. 

The lower solid line if Fig. [4] depicts the the dependence 
of a on / for an ideal gas. The curve deviates from ideal 
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behavior for values of / larger than about 1. At higher 
densities (and pressures) we can improve upon the ideal 
gas behavior by using a virial expansion. From the form 
of the interaction we compute a second virial coefficient 
of B2 = S/tt 2 . As indicated by the upper solid line in 
Fig. 21 inclusion of the second virial coefficient provides 
a good approximation for / up to 3. Clearly there is no 
no simple relation between C and / in this intermediate 
region. 

The focus of this paper was to use the model of hard 
needles to validate the relation between elastic moduli 
and contact probabilities for the exactly solvable model 
of hard needles. However, the model itself has some in- 
teresting features, which will be explored elsewhere [2l| . 
In particular the simplified behavior alluded above in the 



high density limit is related to an incipient critical point. 
The nature and universality of this criticality is related 
to the shapes of the hard objects (in this case, needles). 
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